Field
Embodiments described herein relate generally to a method of reconstructing computed tomography (CT) images, and more specifically to reconstructing CT images using projection data having multiple distinct system matrices geometries.
Description of the Related Art
Computed tomography (CT) systems and methods are widely used, particularly for medical imaging and diagnosis. CT systems generally create images of one or more sectional slices through a subject's body. A radiation source, such as an X-ray source, irradiates the body from one side. A collimator, generally adjacent to the X-ray source, limits the angular extent of the X-ray beam, so that radiation impinging on the body is substantially confined to a planar region (i.e., an X-ray projection plane) defining a cross-sectional slice of the body. At least one detector (and generally many more than one detector) on the opposite side of the body receives radiation transmitted through the body substantially in the projection plane. The attenuation of the radiation that has passed through the body is measured by processing electrical signals received from the detector.
Making projective measurements at a series of different projection angles through the body, a sinogram can be constructed from the projection data, with the spatial dimension of the detector array along one axis (e.g., the horizontal axis) and the time/angle dimension along the other axis (e.g., the vertical axis). The attenuation resulting from a particular volume with the body (e.g., a vertebra) will trace out a sine wave around the spatial dimension along the detector corresponding to the rotation axis of the CT system. Volumes of the body farther from the center of rotation correspond to sine waves with greater amplitudes than those corresponding to volumes nearer the center of rotation. The phase of each sine wave in the sinogram corresponds to the relative angular positions with respect to the rotation axis. Performing an inverse Radon transform (or an equivalent image reconstruction method) reconstructs an image from the projection data in the sinogram, where the reconstructed image corresponding to a cross-sectional slice of the body.
Conventionally, energy-integrating detectors have been used to measure CT projection data. Now, recent technological developments are making photon-counting detectors (PCDs) a feasible alternative to energy-integrating detectors. PCDs have many advantages including their capacity for performing spectral CT. To obtain the spectral nature of the transmitted X-ray data, the PCDs split the X-ray beam into its component energies or spectrum bins and count a number of photons in each of the bins. Since spectral CT involves the detection of transmitted X-rays at two or more energy levels, spectral CT generally includes dual-energy CT by definition.
Many clinical applications can benefit from spectral CT technology, which can provide improvement in material differentiation and beam-hardening correction. Further, semiconductor-based PCDs are a promising candidate for spectral CT, which is capable of providing better spectral information compared with conventional spectral CT technology (e.g., dual-source, kVp-switching, etc.).
One advantage of spectral CT, and spectral X-ray imaging in general, is that materials having atoms with different atomic number Z also have different spectral profiles for attenuation. Thus, by measuring the attenuation at multiple X-ray energies, materials can be distinguished based in the average spectral absorption profile of their constituent atoms (i.e., the effective Z of the material). Distinguishing materials in this manner (i.e., a high-Z material such as bone from a low-Z material such as water) enables a material decomposition from the spectral domain to the material domain. In some instances, this material decomposition is performed using a dual-energy analysis method.
The dual-energy analysis method can be used because the attenuation of X-rays in biological materials is dominated by two physical processes (i.e., photoelectric scattering and Compton scattering). Thus, the attenuation coefficient as a function of energy can be approximated by the decompositionμ(E,x,y)=μPE(E,x,y)+μC(E,x,y),where μPE(E,x,y) is the photoelectric attenuation and μC(E,x,y) is the Compton attenuation. This attenuation coefficient can be rearranged instead into a decomposition of a high-Z material (i.e., material 1) and a low-Z material (i.e., material 2) to becomeμ(E,x,y)≈μ1(E)c1(x,y)+μ2(E)c2(x,y),where c1(x,y) and c2(x,y) are respectively spatial functions describing the concentrations of material 1 and material 2 located at positions (x,y).
Projection data in either the spectral domain or the material domain can be used for image reconstruction. Having obtained and corrected the projection data for measurement and detector artifacts, the projection data is prepared to reconstruct an image of an object OBJ. Generally, the image reconstruction problem can be solved using any of several methods including: back-projection methods (e.g., filtered back-projection), iterative reconstruction methods (e.g., the algebraic reconstruction technique (ART) method and the total variation minimization regularization methods), Fourier-transform-based methods (e.g., direct Fourier method), and statistical methods (e.g., maximum-likelihood expectation-maximization algorithm based methods).
Often the image reconstruction problem will be formulated as a matrix equationAƒ=g,where g are the projection measurements of the X-rays transmitted through an object space that includes the object OBJ, A is the system matrix describing the discretized line integrals (i.e., the Radon transforms) of the X-rays through the object space, and ƒ is the image of object OBJ (i.e., the quantity to be solved for by solving the system matrix equation). The image ƒ is a map of the attenuation as a function of position.
One method of solving the system matrix equation is to recast it as an optimization problem. The optimization problem can also include a regularization term constraining the image solutions to adhere to known image characteristics (e.g., a body will generally have uniform attenuation densities within bodily organs and sharp changes in the attenuation density at the boundary of the organs, thus recommending a total variation minimization regularization term).
When the matrix equation is cast as an optimization problem, any of a number of iterative algorithms can be used to solve the optimization problem. For example, the algebraic reconstruction technique (ART) uses an iterative method of affine projections by refining estimates of the image ƒ through a succession of projections onto each of the affine spaces defined by the rows of the system matrix equation. Performing this series of affine projections and then repeating them again and again causes the image estimate to converge to a solution of the system matrix equation. Interjecting a regularization condition/term between series of affine projections imposes constraints on the solution, and ensures that the reconstructed image converges to a physically meaningful and realistic solution. For example, bodies are known to absorb radiation, making gain through a body (i.e., negative attenuation) an unphysical result. Therefore, a non-negative regularization condition can be imposed to eliminate unphysical gain regions in the reconstructed image.
Iterative reconstruction algorithms augmented with regularization conditions can produce high-quality reconstructions using only a few views even in the presence of significant noise. For few-view, limited-angle, and noisy projection scenarios, the application of regularization operators between reconstruction iterations seeks to tune the final result to be consistent with a certain a priori model. For example, enforcing positivity, as discussed above, is a simple but common regularization scheme.
A second regularization condition is total variation (TV) minimization in conjunction with projection on convex sets (POCS). The TV-minimization algorithm assumes that the image is predominantly uniform over large regions with sharp transitions at the boundaries of the uniform regions. When the a priori model corresponds well to the image object, these regularized iterative reconstruction algorithms can produce impressive images even though the reconstruction problem is significantly underdetermined (e.g., the few-view scenario), missing projection angles, or the data is noisy. The optimization problem with TV regularization can be expressed as
            argmin      u        ⁢          {                                    1            2                    ⁢                                                                  Af                -                g                                                    2            2                          +                  λ          ⁢                                                                  (                                                                        ∇                                                                                  ⁢                    f                                                                    )                                                    1                              }        ,wherein the last term, the l1-norm of the gradient-magnitude image, is the isotropic TV semi-norm. The spatial-vector image ∇u represents a discrete approximation to the image gradient. The expression |∇ƒ| is the gradient-magnitude image, an image array whose pixel values are the gradient magnitude at the pixel location.
Having formulated the optimization problem, there are many choices of algorithms for solving the optimization problem, some of which are more efficient than others, depending on the details of the optimization problem. For example, the dual-primal Chambolle-Pock (CP) algorithm can be efficient for solving certain image reconstruction optimization problems.
Two other related methods for iteratively solving optimization problems are the alternating direction method of multipliers (ADMM) and the augmented Lagrangian multiplier methods. These two methods are examples of splitting-based iterative algorithms. These splitting-based iterative methods exhibit the efficiencies introduced by subdividing a single optimization problem into a series of subproblems that can be solved in an iterative manner.
Although it can be developed in a slightly more general form, the following problem formulation is sufficient for most applications of the ADMM:
            min              x        ∈                  ℝ          n                      ⁢                  ⁢          f      ⁡              (        x        )              +      g    ⁡          (      Mx      )      where, M is an m×n matrix, often assumed to have full column rank, and ƒ and g are convex functions on n and m, respectively. In addition to values in, the functions ƒ and g can also have the value +∞, so that constraints may be “embedded” in ƒ and g, in the sense that if ƒ(x)=∞ or g(Mx)=∞, then the point x is considered to be infeasible.
By appropriate use of infinite values for ƒ or g, a very wide range of convex problems may be modeled by this optimization problem. To make the discussion more concrete, however, herein is described a simple illustrative example that fits very readily into this form without any use of infinite function values, and resembles in basic structure many of applications responsible for the resurgence of interest in the ADMM: the “lasso” or “compressed sensing” problem. This problem takes the form
                    min                  x          ∈                      ℝ            m                              ⁢                        1          2                ⁢                                                        Ax              -              b                                            2                      +          v      ⁢                                  x                          1              ,where A is a p×n matrix, bεP, and ν>0 is a given scalar parameter. The idea of the model is find an approximate solution to the linear equations Ax=b, but with a preference for making the solution vector xεn sparse; the larger the value of the parameter ν, the more the model prefers sparsity of the solution versus accuracy of solving Ax=b. While this model clearly has limitations in terms of finding sparse near-solutions to Ax=b, it serves as a good example application of the ADMM; many other now-popular applications have a similar general form, but may use more complicated norms in place ∥•∥1; for example, in some applications, x is treated as a matrix, and one uses the nuclear norm (the sum of singular values) in the objective to try to induce x to have low rank.
Now, the classical augmented Lagrangian method and the ADMM are described for this optimization problem. First, note that the optimization problem can be rewritten, introducing an additional decision variable zεm asmin ƒ(x)+g(z)ST Mx=z The classical augmented Lagrangian algorithm, is described in the case of the above formulation by the recursions
            (                        x                      k            +            1                          ,                  z                      k            +            1                              )        ∈                            arg          ⁢                                          ⁢          min                                      x            ∈                          ℝ              n                                ,                                          ⁢                      z            ∈                          ℝ              m                                          ⁢              {                              f            ⁡                          (              x              )                                +                      g            ⁡                          (              z              )                                +                      〈                                          λ                k                            ,                              Mx                -                z                                      〉                    +                                                    c                k                            2                        ⁢                                                                            Mx                  -                  z                                                            2                                      }                                ⁢                  λ                  k          +          1                    =                        λ          k                +                                            c              k                        ⁡                          (                                                Mx                                      k                    +                    1                                                  -                                  z                                      k                    +                    1                                                              )                                .                    where, {λk} is a sequence of estimates of the Lagrange multipliers of the constraints Mx=z, while {(xk, zk)} is a sequence of estimates of the solution vectors x and z, and {ck} is a sequence of positive scalar parameters bounded away from 0. The notation a, b denotes the usual Euclidean inner product aTb.
In this setting, the standard augmented Lagrangian algorithm is not very attractive because the minimizations of ƒ and g in the second subproblem in the above expression are strongly coupled through the term
                    c        k            2        ⁢                                    Mx          -          z                            2        ,and hence the subproblems are not likely to be easier to solve than the original problem. In contrast, the alternating direction method of multipliers (ADMM) is formulated into subproblems that are easier to solve. The ADMM problem takes the form
                    ⁢                  x                  k          +          1                    ∈                                    arg            ⁢                                                  ⁢            min                                x            ∈                          ℝ              n                                      ⁢                  {                                    f              ⁡                              (                x                )                                      +                          g              ⁡                              (                                  z                  k                                )                                      +                          〈                                                λ                  k                                ,                                  Mx                  -                                      z                    k                                                              〉                        +                                          c                2                            ⁢                                                                                      Mx                    -                                          z                      k                                                                                        2                                              }                                z              k        +        1              ∈                            arg          ⁢                                          ⁢          min                          z          ∈                      ℝ            m                              ⁢              {                              f            ⁡                          (                              x                                  k                  +                  1                                            )                                +                      g            ⁡                          (              z              )                                +                      〈                                          λ                k                            ,                                                Mx                                      k                    +                    1                                                  -                z                                      〉                    +                                    c              2                        ⁢                                                                                                Mx                                          k                      +                      1                                                        -                  z                                                            2                                      }                                ⁢                  λ                  k          +          1                    =                        λ          k                +                              c            ⁡                          (                                                Mx                                      k                    +                    1                                                  -                                  z                                      k                    +                    1                                                              )                                .                    Clearly, the constant terms g(zk) and ƒ(xk+1) as well as some other constants, may be dropped from the respective minimands in the above expressions. Unlike the classical augmented Lagrangian method, the ADMM essentially decouples the functions ƒ and g. In many situations, this decoupling makes it possible to exploit the individual structure of the ƒ and g so that the above expressions may be computed in an efficient and perhaps highly parallel manner.
The discussion thus far has centered on reconstructing a CT image for a single CT scanner system. A CT scanner system having multiple measurement methods/geometries can be treated identically to the single measurement method/geometry CT system by treating each measurement method and measurement geometry as a separate CT scanner system and reconstructing a separate image for each measurement method and measurement geometry. However, reconstructing separate images ignores the shared information among the multiple measurement methods/geometries that can be used to improve the image quality. Therefore, a better image reconstruction method will collectively process the projection data obtained using the multiple measurement methods/geometries to reconstruct a combined image.